2. Sistemas para el procesamiento de señales¶

2.1. Definición de sistema¶

Hasta ahora hemos realizado análisis de señales, es decir el estudio de las señales y sus propiedades en el dominio del tiempo y frecuencia

En esta unidad nos enfocaremos en el procesamiento de señales, es decir el diseño de sistemas que procesan una señal de entrada y producen una señal de salida

Usaremos

\(x[n]\) para denotar la señal (discreta) de entrada y \(X[k]\) su espectro
\(y[n]\) para denotar la señal (discreta) de salida e \(Y[k]\) su espectro

2.1.1. Ejemplos de sistemas¶

Utilizando sistemas podemos modificar una señal para mejorar su calidad o remover efectos indeseados

Un sistema para reducir el ruido de una señal de electroencefalograma (EEG)

Un sistema para mejorar una imagen fuera de foco (sharpening)

Un sistema para eliminar el eco de un audio

En lo que sigue realizaremos una clasificación de los sistemas según sus propiedades

2.2. Sistemas sin memoria¶

Diremos que un sistema \(\Phi\) es un sistema sin memoria si

\[ y[n] = \Phi(x[n]), \]

es decir que la salida del sistema en un instante \(n\) dado depende solo de la entrada en ese mismo instante

Veamos algunos ejemplos

2.2.1. Sistema atenuador/amplificador ideal¶

\[ y[n] = A x[n], \]

donde \(A>0\) se llama ganancia

Este sistema puede atenuar la entrada si \(0<A<1\) y amplificar si \(A>1\)

2.2.2. Sistema saturador (clamp)¶

\[\begin{split} y[n] = \begin{cases} B &x[n] > B \\x[n] & x[n] \in [-B, B]\\ -B & x[n] < -B\end{cases} \end{split}\]

Este sistema limita los valores de la señal de entrada en un rango fijo

2.2.3. Sistema rectificador¶

\[ y[n] = | x[n] | \]

Este sistema eliminar la parte negativa de la señal de entrada

2.3. Sistema Lineal¶

Diremos que un sistema \(\Phi\) es lineal si cumple con las siguientes propiedades

2.3.1. Homogeneidad¶

Un cambio en la amplitud de la entrada produce un cambio equivalente en la salida

\[ \Phi(cx[n]) = c \Phi(x[n]) = c y[n] \]

2.3.2. Aditividad¶

Señales que se suman en la entrada producen señales que se suman en la salida

\[ \Phi(x_1[n] + x_2[n]) = \Phi(x_1[n]) + \Phi(x_2[n]) = y_1[n] + y_2[n] \]

Es decir que las señales pasan por el sistema sin interactuar entre ellas

2.3.3. Otras propiedades de los sistemas lineales¶

Producto de las propiedades anteriores se tiene que una cascada de sistemas lineales forma un sistema lineal equivalente

Y la cascada de sistemas es conmutativa: El orden de los sistemas en la cascada no altera el resultado final

Los sistemas lineales también cumplen el Principio de superposición

Si descomponemos una señal en \(M\) componentes: \(x[n] = x_1[n] + x_2[n] + \ldots + x_M[n]\)
Y aplicamos un sistema lineal a cada componente: \(y_j[n] = \Phi(x_j[n])\)
Podemos recuperar la salida total usando aditividad: \(y_1[n] + y_2[n] + \ldots + y_M[n] = y[n]\)

2.4. Sistemas con memoria¶

Un sistema \(\Phi\) es un sistema con memoria si su salida actual depende sólo de la entrada actual, las entradas anteriores o las salidas anteriores

\[\begin{split} \begin{align} y[n] = \Phi(x[n], & x[n-1], x[n-2], \ldots, x[0], \\ \nonumber & y[n-1], y[n-2], \ldots, y[0]) \nonumber \end{align} \end{split}\]

esto también se conoce como sistema causal

Un sistema con memoria no-causal usa entradas futuras (es decir \(x[n+1]\), \(x[n+2]\), …) y por ende solo se puede implementar de forma offline, es decir una vez que sea ha observado toda la señal

Veamos algunos ejemplos de sistemas con memoria

2.4.1. Sistema con un retardo (delay)¶

Definido como

\[ y[n] = x[n-m], \]

donde

la salida depende solo de «una» entrada anterior
el valor de m define que tan «antigua» es la entrada pasada

El delay o retarno no afecta la amplitud de los componentes frecuenciales de la señal pero si su fase, como muestra la siguiente animación

%%capture

fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m)) 
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))

def update(m):
    ax[0].cla(); ax[0].plot(n, x(m));
    X = sfft.fftshift(sfft.fft(x(m)))
    Xm = np.absolute(X); Xp = np.angle(X)
    # Espectro de magnitud:
    ax[1].cla(); ax[1].plot(f, Xm); 
    # Espectro de fase enmascarado con el espectro de magnitud
    ax[2].cla(); ax[2].plot(f, Xm*Xp/np.amax(Xm)); 
    ax[2].set_ylim([-np.pi, np.pi])
    angle_delay = Xp[np.argmax(Xm)]
    ax[2].set_title("%0.4f [rad], %0.0f [deg]" % (angle_delay, 180*angle_delay/np.pi))
    return ()

anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=True)

HTML(anim.to_html5_video())

2.4.2. Sistema reverberador o eco¶

Definido como

\[ y[n] = x[n] + A x[n-m], \]

donde

la salida depende de una entrada «pasada» y la entrada actual
la ganancia controla si el «eco» se atenua o amplifica

Al contrario del sistema anterior, el eco si puede modificar el espectro de amplitud.

Notemos el efecto de interferencia constructiva y destructiva al modificar el retardo, como muestra la siguiente animación

%%capture
fig, ax = plt.subplots(3, figsize=(6, 6), tight_layout=True)
n = np.arange(0, 200, step=1)
x = lambda m: np.sin(2.0*np.pi*0.05*(n-m)) 
f = sfft.fftshift(sfft.fftfreq(d=1, n=len(n)))
A = 1.
def update(m):
    y = x(0) + A*x(m)
    X = sfft.fftshift(sfft.fft(y))
    for ax_ in ax:
        ax_.cla()
    ax[0].plot(n, x(0), label='x[n]')
    ax[0].plot(n, A*x(m), label='A*x[n-m]')
    ax[0].legend()
    ax[1].plot(n, y); 
    ax[1].set_ylim([-A-1.1, A+1.1])
    ax[2].plot(f, np.absolute(X)); 
    ax[2].set_ylim([-20, (1+A)*len(n)/2 + 20])

anim = animation.FuncAnimation(fig, update, frames=40, interval=150, blit=False)

HTML(anim.to_html5_video())

Ejemplo de interferencia destructiva en una onda mecánica: https://www.youtube.com/watch?v=IU8xeJlJ0mk

2.4.3. Sistemas con múltiples ecos¶

Pueden combinarse más retardos para hacer un sistema reverberante más complejo

Por ejemplo

\[ y[n] = x[n] + A_1 x[n-m_1] + A_2 x[n-m_2] + A_3 x[n-m_3] + \ldots, \]

Que se escucha como

Fs = 44100; 
n = np.arange(0, 4, step=1.0/Fs) 
x = lambda m: np.sin(2.0*np.pi*880*(n-m))*np.exp(-(n-m)**2/0.5**2)*np.heaviside(n-m, 0)
y = x(0) + 0.5*x(1.) + 0.25*x(2.) + 0.125*x(3.)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(n, y,  line_width=3)
p1.xaxis[0].axis_label = 'Tiempo [s]'
show(p1)
Audio(y, rate=Fs, normalize=False)

2.5. Sistema FIR¶

Generalizando el ejemplo de sistema lineal reverberante a \(L\) retardos llegamos a

\[\begin{split} \begin{align} y[n] &= h[0] x[n] + h[1] x[n-1] + h[2] x[n-2] + \ldots + h[L] x[n-L] \nonumber \\ &= \sum_{j=0}^{L} h[j] x[n-j] \nonumber \\ &= (h* x)[n] \nonumber \end{align} \end{split}\]

que se puede modelar como una convolución discreta entre \(h\) y \(x\)

Este sistema se conoce como

Sistema FIR (finite impulse response)
Sistema MA (moving average)
Sistema todo-zeros

y es de orden L (posee L+1 coeficientes)

2.5.1. Intepretación como media movil (MA)¶

El sistema FIR es equivalente a una media movil ponderada que se aplica sobre la entrada

Los coeficientes del sistema son los ponderadores

Por ejemplo sea un sistema de 3 coeficientes \(h[0]=a\), \(h[1]=b\) y \(h[2]=c\)

\[\begin{split} \begin{align} y[n] = (h*x)[n] &= \sum_{j=0}^{2} h[j] x[n-j] \nonumber \\ &= a x[n] + b x[n-1] + c x[n-2] \nonumber \end{align} \end{split}\]

donde cada salida se calcula a partir de

\[ \overbrace{x[0], x[1], x[2]}^{y[2]} , x[3], x[4], \ldots \]

\[ x[0], \overbrace{x[1], x[2] , x[3]}^{y[3]}, x[4], \ldots \]

\[ x[0], x[1], \overbrace{x[2] , x[3], x[4]}^{y[4]}, \ldots \]

Para obtener el valor de \(y[0]\) e \(y[1]\) se deben establecer «condiciones de borde», como por ejemplo \(x[-2] = x[-1]= 0\)

2.5.2. Ejemplo: Eliminando ruido blanco aditivo¶

Si tenemos una señal contaminada con ruido blanco aditivo podemos usar un sistema FIR promediador para «suavizar la contaminación»

La animación muestra un sistema FIR con \(L\) coeficientes idénticos e iguales a \(1/L\) que se convoluciona con la señal contaminada intentando recuperar la «señal limpia»

%%capture
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C) 
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
# Sistema
L = 10; h = np.ones(shape=(L,))/L; 
# Acumulator
y = np.zeros_like(n, dtype=np.float)

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
ax[0].plot(n, x, 'k.'); 
ylims = ax[0].get_ylim()

def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()

def update(m):
    ax[1].cla(); ax[2].cla()
    c = np.zeros_like(n, dtype=np.float); 
    c[m:m+L] = h
    ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+L])
    ax[2].plot(n, y, label='señal recuperada');  
    ax[2].plot(n, x_clean, 'g', label='señal limpia'); 
    ax[2].legend()
    ax[2].set_ylim(ylims)
    ax[2].plot([m, m], [ylims[0], ylims[1]], 'r--', alpha=0.5)
    return ()

anim = animation.FuncAnimation(fig, update,init_func=init, 
                               frames=100-len(h), interval=150, blit=True)

HTML(anim.to_html5_video())

2.5.3. Ejemplo: Encontrando cambios en una señal¶

Podemos usar un sistema «diferenciador» para detectar cuando una señal cambia de valor como se muestra en la siguiente animación

%%capture

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
n = np.arange(0, 100, step=1)
x = np.zeros_like(n, dtype=np.float)
x[20:] += 1.; x[40:] += 1.; x[80:] -= 1.;
ax[0].plot(n, x)
# System:
h = np.array([-0.5, 0.5])
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()

def update(m):
    c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
    ax[1].cla(); ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+len(h)])
    ax[2].cla(); ax[2].plot(n, y);  
    ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
    return ()

anim = animation.FuncAnimation(fig, update, init_func=init, 
                               frames=100-len(h), interval=150, blit=True)

HTML(anim.to_html5_video())

2.5.4. Ejemplo: Removiendo una tendencia¶

En el ejemplo siguiente tenemos una señal sinusoidal que nos interesa separar de una señal más lenta/suave denominada tendencia

Podemos usar un sistema FIR para eliminar la tendencia

%%capture

fig, ax = plt.subplots(3, figsize=(7, 6), tight_layout=True)
np.random.seed(0); 
n = np.arange(0, 100, step=1)
C = np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/30**2)
x = np.sin(2.0*np.pi*0.1*n) + 5*np.random.multivariate_normal(np.zeros_like(n), C)
ax[0].plot(n, x); ylims = ax[0].get_ylim()
# System:
L = 5; h = -np.ones(shape=(L,))/L; h[L//2] += 1
# Acumulator
y = np.zeros_like(n, dtype=np.float)
def init():
    global y
    y = np.zeros_like(n, dtype=np.float)
    return ()
    
def update(m):
    c = np.zeros_like(n, dtype=np.float); c[m:m+len(h)] = h
    ax[1].cla(); ax[1].plot(n, c); 
    y[m] = np.sum(h*x[m:m+len(h)])
    ax[2].cla(); ax[2].plot(n, y);  
    ax[2].plot([m, m], [-0.5, 0.5], 'r--', alpha=0.5)
    return ()
    
anim = animation.FuncAnimation(fig, update, init_func=init, 
                               frames=100-len(h), interval=200, blit=True)

HTML(anim.to_html5_video())

2.6. Convolución con scipy¶

En la práctica, si queremos convolucionar una señal con nuestro sistema, podemos usar la función de scipy

scipy.signal.convolve(in1, # Señal de entrada
                      in2, # Coeficientes del sistema
                      mode='full', 
                      method='auto'  
                     )

donde el argumento method puede ser

direct: Realiza la convolución en el dominio del tiempo
fft: Realiza la convolución multiplicando los espectros
auto: Se decide automaticamente en base al largo de las señales

y el argumento mode indica donde se hace la convolución

Sea una señal \(x=[a,b,c]\) y un sistema \(h=[d, e]\)

Si uso mode=valid el resultado será \(y=[ad+be, bd+ce]\)
Si uso mode=same el resultado será \(y=[ae, ad+be, bd+ce]\), es decir se agregan ceros al principio de \(x\) tal que \(y\) sea del mismo largo que \(x\)
Si uso mode=full el resultado será \(y=[ae, ad+be, bd+ce, cd]\), es decir se agregan ceros al principio y al final de \(x\)

2.6.1. Ejemplo: Eliminando ruido v 2.0¶

Esta vez usamos scipy.signal.convolve y probamos dos sistemas FIR

coeficientes idénticos e iguales a \(1/L\) (rectangular)
coeficientes que decaen suavemente a cero (hamming)

from bokeh.palettes import Dark2_5 as palette
np.random.seed(0)
n = np.arange(0, 100, step=1)
C = 5*np.exp(-0.5*(n[:, np.newaxis] - n[:, np.newaxis].T)**2/10**2)
# Señal de entrada original
x_clean = np.random.multivariate_normal(np.zeros_like(n), C) 
# Señal de entrada contaminada
x = x_clean + 2*np.random.randn(len(n))
L = 20; 

p = [Figure(plot_width=600, plot_height=230, toolbar_location="below") for k in range(3)]

p[0].scatter(n, x, color='black')
p[2].line(n, x_clean, line_width=2, legend_label = 'señal limpia'); 
for name, h, color in zip(['rect', 'hamming'],
                          [np.ones(shape=(L,)), scipy.signal.hamming(L)],
                          palette):

    h = h/np.sum(h)
    p[1].line(np.arange(L), h, color=color, legend_label=name, line_width=2)
    p[2].line(n, scipy.signal.convolve(x, h, mode='same', method='auto'), 
              color=color, legend_label=name, line_width=2); 
p[2].legend.location = "top_left"
show(column(p))

INFO183 Análisis de sistemas lineales